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Abstract:  numerical  solution  procedure  capable  of  a  high  degree  of 

automation  is  presented  for  calculating  the  flow  field 
associated  with  viscid- inviscid  interaction  on  a  body  of 
revolution.  The  method  makes  use  of  ideas  from 
asymptotic  triple  deck  analysis  of  the  trailing  edge 
problem.  A  simpler  problem  is  solved  in  which  two  layers 
are  postulated,  an  outer  inviscid,  rotational  layer  (^frozen 
vorticity\)  which  largely  governs  the  pressure  in  the 
trailing  edge  region,  and  an  inner  turbulent  layer  which  is 
modeled  in  a  simple  manner.  These  two  layers  are  patched  at 
a  fixed  distance  from  the  wall  and  wake  centerline  rather 
than  asymptotically  matched.  The  method  allows  the  viscid- 
inviscid  interaction  solution  to  be  computed  in  four 
iterations,  none  of  which  require  any  user  intervention.  A 
high  degree  of  automation  is  achieved  in  the  computer  code 
because  of  three  factors:  (1)  the  solution  cycle  is  very 
stable,  (2)  all  calculations  use  a  common  streamwise  grid, 
and  (3)  information  is  transferred  to  successive  stages  of 
the  calculation  by  computer  generated  files.  Numerical 
solutions  are  presented  for  three  bodies  and  the  results 
compared  with  published  experimental  data. 
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Introduction 


The  standard  approach  for  solving  the  axisymmetric  viscid-inviscid 
interaction  problem  is  the  displacement  body  method  of  which  numerous 
examples  have  appeared  in  the  literature  [1-3].  In  this  procedure  the 
boundary- layer  displacement  thickness  is  added  to  the  body  to  obtain  an 
approximation  to  the  strong  interaction  pressure  distribution.  Sines  the 
displacement  thickness  depends  on  knowing  the  pressure  distribution, 
iteration  is  required  to  find  the  solution.  In  the  direct  mode  (pressure 
distribution  given)  the  boundary- layer  solution  invariably  produces  a 
displacement  thickness  which  has  numerical  noise  in  the  strong 
interaction  region.  This  "noise"  is  then  amplified  into  pressure  wiggles 
in  the  potential  flow  solution.  To  obtain  convergence,  numerical 
smoothing  of  the  pressure  distribution  is  required  between  each 
iteration.  This  process  requires  intervention  by  the  user  and  is  thus 
difficult  to  automate.  A  variation  of  the  displacement  body  method  is  to 
solve  the  boundary- layer  portion  in  the  inverse  mode  which  results  in  a 
numerical  smoothing  effect.  An  example  of  this  approach  is  the  work  of 
Carter  and  Womom  [4], 

A  promising  alternative  approach  for  the  solution  of  the 
axisymmetric  viscid-inviscid  interaction  problem  is  to  adopt  the  main 
ideas  of  asymptotic  analysis  as  applied  to  the  two-dimensional  trailing 
edge  problem.  This  approach  has  a  firm  foundation  mathematically  and,  as 
it  turns  out,  can  be  readily  automated. 

The  turbulent  flow  near  a  2-D  trailing  edge  was  first  analyzed  using 
asymptotic  methods  by  Melnick,  Chow  and  Mead  [5]  who  found  that  a  three- 
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layer  structure  develops  which  extends  over  a  streamwise  distance  of  the 
order  of  (1/ln  Re).  In  the  outer  layer,  because  of  the  small  streamwise 
length  scale,  viscous  terms  can  be  neglected  leaving  only  the  inertial 
terms.  Thus  the  flow  in  this  layer  is  inviscid  but  rotational,  the 
vorticity  arising  from  the  upstream  boundary  layer.  An  inner  viscous 
layer  is  required  to  satisfy  the  no-slip  condition  and  an  intermediate, 
or  blending,  layer  accounts  for  the  different  rates  of  growth  of  the 
inner  and  outer  layers.  Although  the  extent  of  the  streamwise  effect  ma; 
differ,  the  same  type  of  three-layer  structure  arises  in  the  vicinity  of 
an  axisymmetric  trailing  edge. 

Analytical  solutions  of  the  2-D  turbulent  trailing  edge  using  the 
three- layer  asymptotic  approach  have  been  obtained  by  Melnick,  Chow  and 
Mead  for  a  cusped  trailing  edge  and  by  Melnick  and  Grossman  [6]  for  a 
small  but  nonzero  trailing  edge  angle.  These  solutions  clearly 
demonstrate  the  importance  of  normal  pressure  gradient  (via  the  outer 
rotational  inviscid  layer)  in  the  correct  description  of  the  flow  at  the 
trailing  edge. 

Prior  to  the  work  of  Melnick  et  al.,  the  only  other  analytical 
treatment  of  the  2-D  turbulent  trailing  edge  was  by  Kuchemann  [7].  He 
considered  a  wedge-shaped  trailing  edge  of  finite  angle  in  the  presence 
of  rotational,  inviscid  flow  with  constant  vorticity.  Since  the  inner 
layers  were  neglected,  a  slip  velocity  at  the  wall  was  required  to  avoid 
a  constant  pressure  on  the  wedge  surface.  His  main  finding  was  that  the 
velocity  field  induced  by  the  vorticity  in  the  wake  retarded  the  flow  on 
the  wedge  upstream  of  the  trailing  edge  which  led  to  separation  if  the 
wedge  angle  was  large  enough. 
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In  the  axisymmetric  case,  Geller  [8]  and  Hoffman  [9]  have  made  use 
of  an  outer  rotational,  inviscid  layer  to  determine  flow  near  the 
trailing  edge.  In  Geller' s  procedure  a  local  body  oriented  Cartesian 
coordinate  system  is  used  and  3v/3x  is  neglected  in  the  vorticity.  As  a 
result,  the  velocity  profile  at  a  given  body  station  is  obtained  by 
solving  an  ordinary  differential  equation  rather  than  an  elliptic  partial 
differential  equation,  as  is  ordinarily  the  case.  He  patches  the 
inviscid,  rotational  profile  at  a  given  normal  distance  from  the  body 
surface  with  a  one-seventh  power  law  turbulent  boundary- layer 
representation  in  the  inner  region.  Agreement  of  calculated  velocity 
profiles  near  the  trailing  edge  with  experimental  data  is  very  good.  He 
also  finds  that  the  computed  velocity  profile  in  the  outer  region  is  not 
sensitive  to  the  assumed  profile  in  the  inner  region.  Hoffman  solves  the 
elliptic  problem  in  terms  of  the  stream  function.  His  solution,  which 
neglects  patching  with  the  inner  layer,  allows  for  velocity  slip  at  the 
wall,  as  in  the  2-D  treatment  by  Kuchemann.  As  a  result  of  this  neglect, 
he  finds  the  solution  to  be  fairly  sensitive  to  choice  of  the  initial 
station  where  the  rotational  calculation  begins. 

In  the  present  work  the  findings  of  2-D  asymptotic  analysis  are 
applied  in  a  approximate  fashion  to  the  axisymmetric.  trailing  edge. 
Instead  of  a  full  triple  deck  treatment,  only  the  outer  layer  is  solved 
exactly,  the  intermediate  layer  is  neglected  and  the  inner  layer  is 
modeled  in  a  simple  manner.  This  work,  therefore,  represents  an 
extension  of  Reference  9. 

Since  asymptotic  matching  is  inappropriate  in  the  present  scheme, 
the  inner  and  outer  layers  are  "patched"  at  a  small,  fixed  distance  from 
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the  body  surface  and  wake  centerline.  The  inner  layer  is  modeled  using  a 
power-law  velocity  profile  to  represent  the  lower  portion  of  the 
turbulent  boundary  layer,  the  same  as  used  by  Geller.  Because  the  wake 
cannot  be  neglected  in  the  elliptic  problem,  a  simplifying  assumption  is 
required  for  its  inner  layer.  The  assumption  is  made  that  close  to  the 
wake  centerline  the  vorticity  varies  linearly  with  radial  distance.  The 
presence  of  the  inner  layer  results  in  a  boundary  condition  of  the  third 
kind  at  the  lower  edge  of  the  outer  layer  (the  "patch  boundary").  Once 
the  solution  in  the  inviscid,  rotational  domain  is  known,  the  pressure  on 
the  body  and  wake  centerline  is  obtained  by  integrating  the  normal 
momentum  equation  downward  from  the  patch  boundary. 

The  reasons  why  the  "frozen  vorticity"  approach  can  be  highly 
automated  and  thus  yield  an  easy-to-use  computer  code  are  (1)  a 
considerable  commonality  exists  in  computational  grids  for  the  inviscid 
irrotational  and  rotational  flowfieids  as  well  as  the  boundary  layer,  (2) 
the  same  body  curve  fit  is  used  for  all  three  flows,  and  (3)  the  starting 
profile  for  the  inviscid  rotational  solution  is  generated  automatically 
using  computer-generated  data  files.  The  boundary- layer  solution  with 
viscid-inviscid  interaction  effects  included  takes  only  four  steps:  (1) 
a  potential  flow  solution,  (2)  a  corresponding  boundary- layer  solution, 
(3)  a  frozen  vorticity  solution,  and  (4)  a  final  boundary- layer  solution. 
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Equations  of  Motion 

In  cylindrical  coordinates,  (x,r,9)  were  3/39  =  0  (axial  symmetry) 
inviscid  potential  or  rotational  flow  is  governed  by  the  equation 


xx 


+  Grr  - 


7Gr  = 


-  r-F ( tp ) 


(1) 


where  a  lower  case  subscript  denotes  partial  differentiation  and  G  is  the 
perturbation  stream  function  defined  by 


G  *  ip  -  ^  r2  . 


(2) 


The  perturbation  stream  function  is  'used  because  it  is  bounded  as  r  ->  ». 
The  reduced  vorticity  F  is  related  to  the  vorticity  C  by 


C  “  rF  ( ip )  . 


(3) 


The  quantity  F  is  uniquely  determined  by  the  upstream  boundary  condition. 
For  potential  flow  F  =  0. 
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where 

A  * 

37  - 

Xr 

(J^)n  +  t1!  . 

(6) 

B  = 

s 

y  ’ 

(7) 

C  * 

a 

Y  * 

(8) 

D  = 

iu=)5 

xn 

-  <JS)n  +  — ]  , 

(9) 

E  = 

X" 

(10) 

and 

J  = 

Jacobian  = 

<xCrn  *  *nr5rl  • 

(11) 

a  = 

2  .  2 

x  +  r 
n  n 

9 

(12) 

0  = 

x5xn  +  r£rn 

9 

(13) 

y  = 

2  2 

x.  +  r.  . 

r  r 

(14) 

•  -  '  k.  _ 

4  t. 


Body  Fitted  Coordinates 

A  C-grid  is  the  appropriate  coordinate  system  in  the  present  problem 
for  proper  flow  resolution  at  the  nose  and  downstream  of  the  tail.  A 
suitable  C-grid,  consisting  of  three  regions,  can  be  generated 
analytically,  as  illustrated  in  Fig.  1.  The  transformations  are  as 
follows: 

(I)  Orthogonal  wrap-around  grid,  0  S  £  S  . 


xb(£)  -  q  sin  <p(£) 

(15) 

rb(£)  +  n  cos  $(£) 

(16) 

where  xb,  rb  and  4>  as  well  as  £  and  q  are  defined  in  Fig.  2. 


Gilbert  H.  Hoffman  Page 

(II)  Sheared  grid,  S  ;  i  ^  . 


dx  =  cos  «(5)dS 

(17 

r  =  rb(0  +  q  . 

(18 

Note  that  for  Regions  I  and  II  to  be  compatible, 

at  the  maximum  body  diameter. 

the  junction  must  occur 

(Ill)  Cartesian  grid,  <  £  <  5^  . 

x  =  5 

(19 

r  =  r|  . 

(20 

The  metric  coefficients  x^,  x^,  etc.  are  obtained  analytically  from 
Eqs.  (15)  through  (20). 

Potential  Flow  Boundary  Cond itions 


For  the  potential  flow  case  Dirichlet  boundary  conditions  are 
specified  on  three  sides  of  the  computational  domain  (consisting  of 
Regions  I,  II  and  III)  and  a  Neumann  condition  on  the  fourth  side.  On 
the  stagnation  line  (left  boundary)  -^  =  0  and  hence 


G(0,n)=0  ,  OSnSriu  . 


(21 
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Since  the  body  and  wake  are  part  of  the  same  streamline  as  the  stagnation 
line,  the  condition  V  =  0  holds  which  leads  to 


GU.O) 


.  1 


-  \  r?  ,  C  i  5  S  . 

ZD  1 


*r  s  €  s  5r 


(22a) 


(22b) 


At  the  outflow  plane  (right  boundary)  the  flow  is  assumed  to  be  uniform 
and  parallel,  the  same  as  the  free  stream,  and  hence 


G(5r,h)  -  0  ,  0  £  n  -  Hu 


Finally,  the  outer  boundary  is  positioned  far  enough  from  the  body  so 
that  the  free-stream  condition  on  velocity  prevails,  viz. 


GnU.Mu)  =  0  ,  0  S  £  S 


Equation  (24)  produces  a  better  behaved  numerical  solution  procedure  that 
its  Dirichlet  counterpart  G(5,nu)  =  0. 


Frozen  Vorticitv  Boundary  Conditions 
In  the  frozen  vorticitv  case  the  computational  domain  consists  cf 
part  of  Region  II  and  all  of  Region  III  with  the  lower  boundary  displaced 
upward  a  distance  A  from  the  body  and  wake  centerline.  The  boundary 
condition  on  the  upper  boundary  n  =  n^j  is  the  same  as  in  the  potential 
flow  case.  On  the  left  or  inflow  boundary  the  stream  function  and 
vorticity  are  specified  as  a  composite  function  of  the  potential  flow  and 
boundary- layer  solutions,  viz. 


*  -  aw  n  m  w  aw  xw  a^.  «  »  w*  v  w  w  j  trm  r  _  <-w  -y-_  ^  *4"^  >r  i »-«*->  a .  “»■  ^j*  w^f  _w  *-^r  "ji  — ^  ^ 


Gilbert  H.  Hoffman 


G(5j-\/,n)  3  Gc(n) 


FUpy.n)  *  Fc(n) 


Page  10 


0  £  n  £  no  . 


On  the  outflow  plane  (right  boundary)  parallel  flow  is  assumed,  as  in  the 
potential  flow  case.  This  assumption  leads  to 


cp^R’"^  =  ®  » 


plus  the  following  quadrature  relation  between  r  and  p  : 


If  i 

tp  =  -  dr-  on  5R  ,  0  ^  n  -  Hu  . 


where  from  Bernoulli's  equation  and  Eq.  (26) 


u  =  yJcp0lvJ 


The  numerical  details  of  evaluating  the  inflow  and  outflow  boundary 
conditions  will  be  given  later. 

The  boundary  condition  on  the  lower  ("patch")  boundary  is  obtained 


from  considerations  in  the  modeled  turbulent  boundary  layer.  The  assumed 
"power  law"  velocity  profile  in  this  layer  may  be  written  as 


,  v.  n 

u  =  u,  (7) 


where  u,  is  the  axial  velocity  at  v  =  A  and  is  determined  from  the  frozen 


vorticitv  solution.  The  stream  function  and  axial  velocity  are  related 


\  s 


TanTTT 
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Integration  of  Eq.  (30),  using  Eq.  (29),  between  y  =  0  and  v  =  d  yields 


=  vA’(m;rb  +  r73; A)  •  (: 

With  the  aid  of  Eqs.  (2)  and  (30),  is  found  to  be  related  to  the  n- 

derivative  of  G  by 

(G.)i 

U.  =  1  T  -  . 

-1  r , 


The  required  boundary  condition,  upstream  of  the  tail,  is  obtained  by 
combining  Eq.  (2),  evaluated  at  y  =  A  ,  with  Eqs.  (31)  and  (32).  The 
result  is 


GA  '  |  =  *  -  \  rA2  ,  Spy  S  5  S  5T 

A 


where 


K  *  a'(  rhr  rb +  rrs a)  - 


rA  =  rb  +  A  • 

In  the  wake  the  assumption  is  made  that  near  y  =  0  the  vorticity  is 


linear  in  y. 


C  "  ^  ’  I 


Using  Eq.  (36)  and  neglecting  3v/ 3x  (in  keeping  with  the  boundary- layer 
approximation),  the  vorticity  definition  can  be  integrated  to  yield  the 
following  expression  for  the  centerline  velocity: 


=  ud  +  2 


4 JIM ■ Jg J* ji  V 
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The  integral  of  Eq.  (2)  between  y  *  0  and  y  =  A  ,  making  use  of  Has .  (3) 
and  (37),  gives 

|  -t-  |  F(ti/)-d4  .  (33) 

When  rewritten  in  terms  of  and  (G^)^  the  wake  boundary  condition  on 
y  =  4  is  obtained,  viz. 


E  >  E- 


Numerical  Algorithm 


The  main  features  of  the  numerical  method  used  to  solve  the 
potential  flow/frozen  vorticity  equation  are  as  follows: 

(1)  The  transformed  vorticity  equation,  Eq.  (5),  is  written  as  a  first- 
order  system. 

(2)  A  fourth-order  accurate  spline,  sHl.O),  is  used  in  the  n  direction 
to  resolve  the  vortical  layer  with  as  few  nodal  points  as  possible. 

(3)  Second-order  accurate  finite  difference  formulas  are  used  in  the 
^-direction. 

(4)  A  nonuniform  grid  in  n  is  generated  by  the  use  of  a  stretching 
function. 

(5)  The  resulting  system  of  algebraic  equations  is  solved  by  SLCR, 
sweeping  in  the  E-direction. 

To  write  Eq.  (5)  as  a  first-order  system  the  following  auxiliary 


variable  is  defined: 
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H^  -  AH  -  2BH5  +  CG^r  -t-  DGr  =*  -  E 


(41 


The  next  step  is  the  definition  of  the  following  spline  first 
derivatives : 


Then  the  governing  equations,  (40) 

lG  -  H  =  0  , 


(42 

(43 

and  (41),  become 

(44 


and 


-  AH  -  23Hr  +  CGrr  +  DGr  +  E  =  0 


(45 


The  finite  difference  expressions  used  for  the  ^-derivatives  are  those 
given  by  Blottner  [10]  for  a  nonuniform  grid,  namely 


Gi+l,j  ~  Gi-l,j 
5i+l  *  H-l 


and 


j,»  i 
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L»J  5i+l  ‘  5i-l 


’i+l.j 


Gi.j  '  Gi-l,j 
5i  -  Si-i  J 


plus  an  exDression  similar  to  Ea.  (46)  for  ( He- ) ^  4.  The  5  discretizati< 

S  i  >  J 

of  Eq.  (45)  leads  to  the  following  expression  for  unknowns  at  point 


(i,j): 


a.i  4G4  4  -  4H4  4  =  R4  4 

4-»J  -»J  1>J  *»J 


^i+1  '  5i-i  5i+l  ‘  S 


i  +  «i  -  Si-i)  Ci»j 


bi,j  =  Ai,j 


Ri’j  "  €i+i  -  «i-i  K  «i  -  5i-i  '  Di’j}  Gi-L*j 


5i+i  -  5i  +Di*j}  Gi+l-j]  +  5i+l  -  5i-i 


EU  * 


The  unknowns  at  node  point  (i,j)  are  G,  ZG,  H  and  £**.  Therefore,  to 
complete  the  system  two  spline  relations  are  needed.  The  governing 
equation  has  been  written  in  first-order  form  so  that  the  same  spline 
relation,  SL(4,0),  can  be  used  twice.  The  expression  for  SH4.0)  is 
(with  i  subscript  understood)  [11] 


+  <r2i$j-l  +  BBj*j  +  d  +  a)2Z$j  +  CC4$j+1  +  Z$j+1  =  0,  (52) 


where  <j>  denotes  G  or  H,  and 


2'j~(2  +  a) 

An  4.1  (1  +  a)  * 


(53) 


'mmaiiimyiTv  wvav  tv  wm  »fv  vj  w  ifjvwv.W  ’V**.1  ’A’J  *.»  TW.V  *yny  «VT  c*  it  'XVwV.v^.v.v.v  v  vw  w  w  -  v  v;  ^  v.  yjv:  -m 
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BB< 


CCj 


2(1  +  2a) 
Anj.]_(l  +  a)a 


and 


An. 


a  *  C:  = 


J  Anj.i  ’ 


Anj  =  nj+1  -  Hj  • 
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(54) 

(55) 


(55) 


(57) 


The  number  of  unknowns  at  node  point  (i,j)  is  reduced  from  four  to 
two  by  solving  Eqs.  (44)  and  (45)  for  and  then  substituting 

into  the  two  tridiagonal  spline  relations.  The  resulting  set  may  be 
written  as  the  following  tridiagonal  matrix  equation: 


®i,j2i,j-l  +  Ai,j2i,j  +  ^i,j2i,j+l  =  ^i,j»  2  S  j  S  M  ,  (58) 

where  the  two-component  column  vector  of  unknowns  Z^j  is  defined  by 

zTi,j  -  C=.aii,j  ■  (59) 

The  2x2  matrix  elements  A,  B  and  C  and  the  two-component  column  vector 
D  contain  known  quantities.  These  elements  and  components  can  be 
obtained  from  a  combination  of  Eqs.  (44),  (45)  and  (52). 

The  boundary  condition  on  the  lower  boundary  (n=0orn=A)  can  be 
written  in  the  general  form,  valid  for  both  potential  and  frozen 
vorticity  cases, 


Ll- 
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where  a  and  3  are,  except  in  one  instance,  functions  of  £.  The  exception 
is  Eq.  (39)  where  3  depends  on  G. 

The  boundary  conditions  on  the  lower  and  upper  boundaries  together 
with  a  two  point  spline  relation  at  each  boundary  close  the  system  of 
equations  at  line  i.  The  two-point  spline  relation  used  here  is  Eq.  (16) 
of  Reference  11,  which  in  the  present  variables  is 


An 

—  Hi,2  + 


iH  0  (61) 

1  f  i 


and 


Gi,N+l  ‘  Gi,N 


AOfj  An^j  Anjq  „  ap^  _ 

~2~  Hi,N  '  ~2~  Hi,M+l  +  ~  *i,N+l  ~  IT"  •  li,M  =  0  * 

(62) 


The  above  equations  are  fifth-order  accurate  in  Anj  and  have  been  found 
in  the  present  application  to  maintain  better  overall  solution  accuracy 
then  their  lower-order  cousins. 

Upon  elimination  of  using  Eq.  (45),  the  following  matrix 
relations  result  at  the  boundaries: 


Ai,lzi,l  +  ci,lzi,2  3  Di,  1 


(63) 


and 


Bi,N+lzi,N  +  Ai,N+lzi,N+l  3  Di,N+l  •  (64) 

For  potential  flow  the  singularity  at  the  tail  point  and  on  the  wake 
centerline  requires  special  treatment.  As  can  be  seen,  Eq.  (5)  at  r  =  0 
reduces  to  G^  *  0.  This  condition,  written  as  H^j  a  0,  replaces  spline 


relation  (61). 

Along  any  line  i  =  constant  the  set  of  block  tridiagonal  equations 
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For  potential  flow  the  singularity  at  the  tail  point  and  on  the  wake 
centerline  requires  special  treatment.  As  can  be  seen,  Eq.  (5)  at  r  =  0 
reduces  to  G^  =  0.  This  condition,  written  as  =  0,  replaces  spline 

relation  (61). 

Along  any  line  i  =  constant  the  set  of  block  tridiagonal  equations 
(58),  (63)  and  (64)  is  solved  by  L-U  decomposition  using  case  (i)  of 
Reference  12. 


Mao  Junction  Lines 


For  the  coordinate  system  used  in  the  present  problem,  the  metric 
coefficients  will  be  discontinuous  at  the  junction  of  Regions  I  and  II 
and  the  junction  of  Regions  II  and  III.  When  a  body  has  a  series  of 
curvature  discontinuities,  then  Regions  I  and  II  are  re-zoned  with 
additional  map  junctions  introduced  at  each  curvature  discontinuity.  A 
"mapping  region"  is  defined  here  to  mean  that  portion  of  a  map  described 
by  a  coordinate  transformation,  Eq.  (4),  with  continuous  first  and  second 
derivatives. 

At  map  junctions  where  the  metric  coefficients  are  discontinuous, 

Eq.  (48)  must  be  modified  to  account  for  such  discontinuities.  This 
process  is  performed  using  the  generalized  Chmieiewski-Hof fman  method  of 
Reference  13.  The  author  has  found  that  ignoring  these  discontinuities 
leads  to  errors  as  large  as  37  percent  in  the  potential  flow  pressure 
coefficient. 

In  the  C-H  method,  each  adjoining  domain  is  extended  one  step  into 
the  other  to  form  a  line  of  fictitious  nodes.  The  equation  of  motion. 
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Eq.  (5),  is  then  written  separately  in  the  left  and  right  regions  at  the 
junction.  A  condition  of  smoothness  of  derivatives  of  the  solution  must 
also  exist  across  the  junction.  By  combining  the  left  and  right  forms  of 
the  governing  equation  with  the  smoothness  condition,  a  single  spline- 
finite-difference  (SFD)  equation  is  obtained  at  the  junction  line  which 
accounts  for  the  discontinuity  in  metric  coefficients.  For  the  algebraic 
details,  the  reader  is  referred  to  Reference  9. 


,s 

3 


Inflow  and  Outflow  Conditions 


For  the  frozen  vorticity  case  the  conditions  on  the  inflow  boundary 
are  extremely  important  in  determining  the  behavior  of  the  solution  near 
the  body  and  wake  centerline.  Thus,  care  must  be  taken  to  determine 
these  conditions  accurately. 

As  already  mentioned,  these  conditions  are  a  composite  of  the 
potential  flow  and  boundary- layer  solutions  at  the  initial  line.  This 
composite  is  calculated  as  follows: 

(1)  The  boundary- layer  velocity  and  total  head  profiles  for  A  -  n  -  Hjj 
are  merged  smoothly  with  their  potential  flow  counterparts.  For  the 
velocity,  smooth  merging  occurs  in  a  natural  manner.  For  the  total  head, 
the  merging  point  is  taken  where  cpq  in  the  boundary  layer  just  becomes 
unity.  As  y  increases,  the  total  head  in  the  boundary  layer  will  exceed 
unity  because  the  vertical  component  of  velocity  continues  to  increase. 


W 


.■ 


a 


N.N’.S. 


.% 

i 


| 

A 


Gilbert  H.  Hoffman 


Page  19 


(2)  With  u  known  on  the  inflow  boundary  (denoted  by  u-q-y) ,  the  stream 
function  is  determined  by  numerical  integration  of 


The  integrations  is  carried  out  using  the  trapezoidal  rule  formula: 


#FVj  =  ^FVj  _  i  -r 


( UFV j , i  +  uFVj )  , 


where 


A  =  r2 


dAj  =  Aj  -  Aj  _ 


The  perturbation  stream  function  on  the  initial  line  is  computed  from 
Eq.  (2). 

(3)  With  (cPo)pv  known,  the  reduced  vorticity  distribution,  Fjy(r),  is 
determined  from  the  relation: 


l  3(Cpo>] 


where  the  derivative  of  (cPq)vv  with  respect  to  r  is  calculated  using  a 
three-point  unequally  spaced  finite  difference  formula. 

The  outflow  stream  function  distribution  is  also  calculated  using 
the  trapezoidal  rule,  viz. 


^Rj  -  *Rj_1  + 


(uRj.!  +  uRj ) » 
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where  depends  on  ^  through  Eq.  (28).  Thus  for  each  integration  step, 
iteration  is  required  to  determine  ^Rj  during  which  uRj  is  allowed  to  lag 
one  cycle.  The  foregoing  procedure  is  much  simpler  and  produces  results 
almost  as  accurate  as  solving  for  from  a  two-point  boundary  value 
problem  using  spline  discretization. 


Pressure  on  the  Bodv  Surface  and  wake  Centerline 


For  the  frozen  vorticity  case  the  normal  momentum  equation,  with 
viscous  terms  neglected,  is  used  to  determine  the  pressure  on  the  body 
and  wake  centerline.  Upon  transformation  to  (£,n)  coordinates  and  noting 
that  in  Regions  II  and  III  x  =  x(5)»  this  equation  becomes 


rn  ov 


=  2(v  -  —  u)  —  +  2  —  u  — 

Xc  on  Xc  Or 


The  finite  difference  form  of  Eq.  (70)  is  obtained  by  centering  at  5  = 
and  n  =  A/2  and  approximating  p-derivatives  by  central  differences. 

Then,  applying  the  no-slip  conditions  at  the  wall  and  using  eqs.  (17)  and 
(18)  to  determine  the  metric  coefficients,  the  expression  for  the  wall 
pressure  is  found  to  be 


CPi,W  "  CPi,d  +  °"Pi 


where 


'Pi  “  (v’rbu>i,A  vi , A  +  2  Ui.A  J  1  +  rbl  (vE}i,A  ' 


'b  dx 
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The  derivative  of  v  with  respect  to  5  is  computed  by  a  second-order 
accurate  three-point  formula  for  unevenly  spaced  points.  At  the  ends  of 
the  interval,  second-order  accurate  one-sided  formulas  are  used. 

Results  and  Discussion 

The  main  advantage  of  the  frozen  vorticity  approximation  over  its 
competitors  for  solving  the  viscid-inviscid  interaction  problem  is  the 
high  degree  of  automation  possible.  The  key  to  automation  is 
communication  through  automatically  generated  data  files  between  the  two 
computer  codes  that  perform  the  various  calculations.  The  two  codes  are: 

AXFL04  -  This  code  performs  potential  flow-frozen  vorticity 
calculations  using  the  method  given  in  this  paper. 

BL20  -  This  code  computes  laminar/ turbulent  boundary  layers  using 
the  Keller  Box  Method  [14],  The  algebraic  turbulence  model 
of  Reference  15  is  used. 

The  steps  in  a  complete  viscid-inviscid  interaction  calculation  are 
as  follows: 

(1)  Using  AXFL04  in  the  potential  flow  mode,  a  body  pressure 
distribution  is  generated.  Two  data  files  are  created,  one  for  BL20 
containing  the  pressure  distribution  and  body  curve  fit,  and  the  other 
for  AXFLQ4  (frozen  vorticity  mode)  containing  the  potential  flow  velocity 
profile  at  the  initial  value  line. 
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(2)  BL20  is  run  next  to  obtain  a  boundary- layer  solution  using  the 
potential  flow  pressure  distribution  from  AXFLQ4  which  contains  the 
boundary- layer  velocity  profile  at  the  initial  value  line.  A  double 
interpolation  procedure  is  used  to  obtain  values  in  cylindrical 
coordinates  from  those  in  boundary- layer  coordinates. 

(3)  AXFL04  is  now  run  in  the  frozen  vorticity  mode  to  compute  a  modified 
body  pressure  distribution  in  the  viscid-inviscid  interaction  region. 

The  code  automatically  forms  a  composite  initial  value  profile  using  data 
from  steps  1  and  2.  A  data  file  is  created  for  BL20  which  contains  the 
viscid-inviscid  interaction  body  pressure  distribution. 

(4)  Using  the  pressure  distribution  from  step  3,  BL20  is  re-run  for  the 
viscid-inviscid  interaction  body  pressure  distribution. 

A  nonuniform  point  distribution  in  n  is  generated  using  a  one-sided 
Vinokur  stretching  function  [ 16 ] .  The  stretching  function  of  Vinokur  was 
chosen  because  it  produces  a  grid  with  a  uniform  truncation  error 
independent  of  the  governing  equation  or  difference  algorithm. 

Although  the  same  ^-distribution  is  used  in  the  potential  and  frozen 
vorticity  calculations,  the  n  distributions  differ  considerably.  For 
proper  resolution  in  the  frozen  vorticity  case,  about  a  third  of  the  grid 
points  must  be  placed  in  the  thin  vortical  layer  (the  same  thickness  as 
the  boundary  layer) .  To  do  this  requires  a  rapidly  expanding  grid  in  n 
since  typically  ~  2.  This  is  the  main  reason  splines  rather  than 
finite  differences  were  chosen  to  approximate  flow  derivatives  in  n* 

An  accurate  body  curve  fit,  a  prerequisite  in  obtaining  an  accurate 
potential,  boundary- layer  and  frozen  vorticity  solution,  is  achieved  as 
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follows:  With  the  body  shape  rb(x)  given  either  analytically  or 
discretely,  the  corresponding  arc  length  distribution  ?  is  computed  using 
the  spline  formula 

hi  hi 

5i  *  «i-i  +  “  (“i-i  +  ut)  ‘  xi(mi  ‘  mi-0  »  (73> 


where  the  two  cases  that  oust  be  considered  are  shown  in  the  following 
table: 


i 

1  1  1 

case 

range 

u 

m 

1 

1 

h  i 

xb  -  xb(r) 

|  i 

|rb 

>  1 

(1  +  x’2)l/2 

■ 

Arb 

u 

rb  *  rb(x) 

1  • 

0  - 1  rb 

>  i 

U  «• 

’  ,, 

rbrb 

bxb 

i 

u 

Table  1.  Body  Arc  Length  Parameters 

The  first  case,  x  =  x(rb),  is  appropriate  in  the  nose  region  of  a  blunt- 

1  t  t  !  t  1 

nose  body.  The  quantities  xb,  xb,  rb  and  r-D,  required  in  the 

calculations,  are  computed  in  terms  of  arc  length  derivatives  of  these 
quantities  which  in  turn  are  computed  using  a  three-point  uneven  spacing 


formula. 
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Since  depends  on  Eq.  (73)  must  be  solved  at  each  step  by 
iteration,  convergence  is  very  rapid,  usually  requiring  about  four 
cycles.  The  use  of  Eq.  (73)  will  produce  a  5  distribution  accurate  to 
about  four  decimal  places  which  assures  that  the  derivatives  required  in 
the  mapping  are  accurate  and  smooth.  Test  calculations  have  shown  that 
Eq.  (73)  produces  a  5  distribution  three  to  four  times  more  accurate  than 

I 

the  chord  formula  in  regions  where  r^  is  large  and  changing  rapidly.  In 
the  process  of  determining  £,  the  various  mapping  derivatives  involving 
r^  will  automatically  be  determined.  For  the  boundary- layer  solution, 
the  derivative  dua/d£  is  computed  in  the  same  way  as  body  derivatives, 
using  a  three-point  formula  with  unequal  spacing. 

The  three  bodies  used  as  test  cases  in  this  paper  are  the  same  ones 
used  in  Reference  15,  namely,  the  F-57  low-drag  body  of  Parsons  and 
Goodson  [17J,  the  modified  spheroid  of  Patel  [18]  and  the  NSRDC 
streamlined  body  number  1  [19].  These  bodies  were  chosen  because  of  the 
high  quality  of  the  experimental  data  available  for  each. 

The  grid  parameters  used  in  the  SFD  solutions  presented  here  are 
given  in  Tables  2  and  3  for  potential  flow  and  frozen  vorticity  cases 
respectively.  Reynolds  number  for  each  test  body  (based  on  chord  length) 
appears  in  Table  3.  These  Reynolds  numbers  correspond  to  values  given  in 
References  18-20. 
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1 

F-57 

Spheroid 

NSRDC  No.  1 

XR 

1.5 

1.5 

1.5 

2.0 

2.0 

2.0 

49 

54 

104 

35 

39 

78 

Nn 

20 

20 

20 

An^ 

0.01119 

0.01119 

0.01119 

1 

Table  2.  Potential  Flow  Solution  Parameters 
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! 

F-57 

Spheroid 

i 

NSRDC  No.  1  | 

| 

XFV 

0.707 

0.631 

1 

0.537 

XR 

1.5 

1.5 

1.5 

nU 

2.0 

2.0 

2.0 

Nr 

32 

31 

51  I 

14 

16 

25 

Nh 

40 

40 

40 

Aril 

0.00055 

0.00055 

0.00055 

A 

0.005 

0.005 

0.003 

Re 

1.2  x  106 

1.262  x  106 

6.6  x  106 

Table  3.  Frozen  Vorticity  Solution  Parameters 


The  need  for  splines  in  the  n-direction  is  clearly  seen  by  contrasting 
values  for  Ani  in  Tables  2  and  3.  The  frozen  vorticity  value  is  one 
twentieth  as  large  as  the  potential  flow  value. 

The  number  of  grid  points  used  for  the  NSRDC  body  is  about  twice 
that  used  for  each  of  the  others.  The  reason  for  this  was  to  obtain 
adequate  definition  of  the  body  pressure  in  the  tail  region  where  the 
NSRDC  body  has  a  double  reflex. 

Relaxation  factors  of  1.5  and  1.0  were  used  on  G  and  H  in  all 
calculations.  No  optimization  was  attempted  for  these  factors. 
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Convergence  was  considered  accomplished  in  both  potential  and  frozen 
vorticity  runs  when  the  absolute  maximum  change  in  H  was  less  than 
5  x  10  The  number  of  iterations  required  to  reach  convergence  varied 

from  90  to  110  with  typical  CPU  times  on  a  VAX  11/780  computer  ranging 
from  two  to  three  minutes.  In  all  cases  the  maximum  residual  decreased 
smoothly  as  the  number  of  iterations  increased. 

Boundary- layer  solutions  used  the  same  £  spacing  as  the  potential 
flow  calculations  (as  already  mentioned).  The  number  of  points  in  the 
normal  direction  ranged  from  30  at  the  nose  to  about  50  near  the  tail. 

To  resolve  the  turbulent  boundary  layer  properly,  points  were 
nonuniformly  spaced  using  a  geometric  progress  on,  with  the  smallest  mesh 
at  the  wall.  Typical  CPU  times  were  about  45  seconds.  Thus,  a  complete 
viscid-inviscid  solution  cycle  of  four  runs  (potential  flow,  boundary 
layer,  frozen  vorticity  and  final  boundary  layer)  took  at  most  seven  and 
a  half  minutes  of  CPU. 

The  frozen  vorticity  solution,  having  a  truncated  domain  of 
dependence,  will  depend  on  the  location  of  its  initial  value  line,  i.e., 
on  the  parameter  Xjy.  Experience  has  shown  that  at  Reynolds  numbers 
between  10^  and  10^,  (the  range  covered  here)  that  the  viscid-inviscid 
interaction  region  begins  about  70  percent  of  the  chord  back  from  the 
nose.  Thus  a  suitable  choice  for  Xvy  is  about  0.7.  Since  the  onset  of 
the  interaction  region  is  gradual,  the  solution  should  not  depend 
strongly  on  Xp-y.  To  show  that  this  is  the  case,  different  values  of  xpy 
were  tried  for  two  of  the  three  test  bodies.  The  change  in  the  maximum 
body  pressure  coefficient  (which  occurs  near  the  tail)  was  taken  as  a 
measure  of  the  sensitivity  of  the  solution  to  changes  in  xpy.  Tor  the 
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Modified  Spheroid,  a  change  in  Xpy  from  0.681  to  0.716  produced  a  change 


in  maximum  Cp  of  1.9  percent  while  for  the  NSRDC  body,  a  change  in  xcy 
from  0.687  to  0.747  resulted  in  a  change  of  2.0  percent. 

Since  the  solution  is  relatively  insensitive  to  xj-y,  the  implication 
is  that  the  reduced  vorticity  distribution  on  the  initial  value  surface 
is  essentially  "frozen".  This  is  indeed  the  case  and  is  illustrated  by 
Figure  3.  The  results  in  this  figure  are  for  the  NSRDC  body,  but  are 


b 


typical  for  the  other  cases  studied.  The  figure  shows  the  variation  of 
reduced  vorticity  F  with  stream  function  $  for  two  choices  of  Xjy  (0.687 
and  0.747). 

The  "patching"  of  the  inner  and  outer  layers  is  responsible  for 
introducing  the  parameter  A  into  frozen  vorticity  solution.  But  what  is 
A?  Is  it  merely  a  calibration  parameter  or  does  it  correspond  to  some 
physical  quantity  in  the  turbulent  boundary  layer?  To  try  to  find  the 
answer  to  this  question,  A  was  first  assumed  to  be  a  calibration  constant 
dependent  only  on  Reynolds  number.  The  calibration  was  performed  using 


the  F-57  body  by  requiring  the  predicted  maximum  body  Cp  to  be  the  same 


as  the  experimental  value.  Using  this  procedure  A  was  found  to  be  0.005. 
The  same  value  was  also  used  to  compute  the  solution  for  the  Modified 
Spheroid  (since  the  Reynolds  numbers  are  nearly  the  sane).  The  value  of 
A  was  reduced  to  0.003  for  the  NSRDC  body  to  account  for  the  thinning  of 
the  boundary  layer  with  increasing  Reynolds  number.  At  the  initial  value 
line  the  previously  determined  values  of  A  in  law-of -the-wall  units 
varied  from  200  to  750  whereas  in  terns  of  the  boundary- layer  thickness 
the  range  was  from  0.15  to  0.20.  When  applied  to  the  computed  turbulent 
boundary  layer  at  the  initial  plane,  the  calibrated  values  of  A  were 
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found  not  to  correspond  in  any  consistent  manner  to  a  meaningful  length 
(such  as  the  location  of  maximum  Reynolds  stress  or  the  outer  edge  of  the 
logarithmic  region).  Thus  A  seems  merely  to  serve  as  a  model  calibration 
parameter. 

As  an  indication  of  the  sensitivity  of  the  frozen  vorticity  solution 
to  changes  in  A,  several  values  were  run  for  each  body.  Table  A  is  a 
summary  of  the  maximum  pressure  coefficient  for  these  runs.  The  table 
shows  that  for  two  of  the  three  test  cases  the  solution  is  quite 
sensitive  to  variations  in  A.  As  A  decreases  maximum  Cp  is  seen  to 
increase  with  the  percentage  increase  depending  on  the  type  of  body.  The 
Spheroid  is  by  far  the  least  sensitive. 


Body 

A 

CPmax 

x-Location 

F-57 

0.002 

0.1986 

0.889 

F-57 

0.003 

0.1701 

0.889 

F-57 

0.005 

0.1613 

0.889 

Spheroid 

0.002 

0.2195 

0.978 

Spheroid 

0.003 

0.2184 

0.978 

Spheroid 

0.005 

0.2225 

0.978 

NSRDC 

0.001 

0.1793 

0.968 

NSRDC 

0.003 

0.1606 

0.968 

Table  4.  Sensitivity  of  Computed  cpmax  to  A. 
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The  computed  body  pressure  distributions,  both  potential  and  frozen 
vorticity,  are  compared  with  experiment  in  Figures  4-6  for  the  F-57, 
Modified  Spheroid  and  NSRDC  bodies  respectively.  The  values  of  4  used 
are  those  given  previously  in  this  paper.  In  all  cases  the  frozen 
vorticity  pressure  closely  follows  the  experimental  distribution.  For 
the  Spheroid  the  computed  maximum  Cp  is  too  high.  This  is  also  the  case 
in  Reference  3  which  uses  the  displacement  body  method  (the  present  value 
is  0.223,  Reference  3  gives  0.210  and  experiment,  from  Reference  18  is 
0.180). 

Another  item  of  interest  in  the  present  calculation  scheme  is  the 
magnitude  of  the  pressure  change  across  the  modeled  turbulent  layer. 

This  change  has  been  found  to  be  quite  small.  In  a  typical  run  (the  F-57 
body  with  A  =  0.003)  the  increment,  in  terms  of  the  surface  pressure, 
reached  a  maximum  of  1.8  percent  very  close  to  the  tail. 

The  calculations  of  Geller  [8],  using  a  simplified  frozen  vorticity 
technique,  show  that  the  computed  velocity  profile  (for  the  F-57  body) 
and  the  experimental  profile  develop  almost  identically  with  downstream 
distance.  To  check  this  observation,  velocity  profiles  (frozen  vorticity 
and  turbulent  boundary  layer)  were  computed  with  the  present  procedure  at 
several  downstream  stations.  Since  cylindrical  coordinates  are  used  in 
the  comparison,  a  double  interpolation  procedure  was  necessary  to  obtain 
the  boundary- layer  profiles.  The  comparisons  are  shown  for  the  F-57  in 
Figures  7-9  at  three  successive  stations,  x  =  0.844,  0.933  and  0.978. 

The  frozen  vorticity  profile  is  seen  to  evolve  differently  from  its 
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turbulent  counterpart,  becoming  ever  more  concave  as  x  increases. 
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Apparently  the  poorer  agreement  with  experiment  in  the  present  case  is  a 
result  of  the  full  elliptic  problem  being  solved. 

Concluding  Remarks 

A  method  for  solving  the  axisymmetric  viscid- inviscid  interaction 
problem  has  been  presented  which  predicts  body  pressure  distributions 
with  considerable  accuracy  and  which  is  capable  of  a  high  degree  of 
automation.  In  the  method  the  flowfield  is  divided  into  an  outer 
inviscid,  rotational  region  where  the  vorticity  is  "frozen"  and  an  inner 
viscous/ turbulent  region  which  is  modeled  by  a  power- law  velocity 
profile.  The  solution  to  the  frozen  vorticity  problem  gives  the  viscid- 
inviscid  interaction  body  pressure  distribution  which  is  then  used  as 
input  to  the  standard  axisymmetric  turbulent  boundary- layer  equations. 
With  this  formulation,  a  complete  solution  to  the  viscid-inviscid 
interaction  problem  is  obtained  in  four  computation  cycles.  The  method 
requires  no  user  intervention  between  cycles  and  is  more  economical  than 
solving  the  complete  Navier-Stokes  equations  or  the  displacement  body 
problem. 

Results  with  the  present  method  show  that  the  outer  rotational  layer 
(frozen  vorticity  solution)  controls  the  pressure  in  the  viscid-inviscid 
interaction  region  near  a  body  surface.  This  observation  is  in  agreement 
with  the  findings  of  triple  deck  theory. 

The  present  method  is  an  approximate  application  of  a  triple  deck  in 
which  the  middle  layer  is  ignored  and,  instead  of  asymptotic  matching, 
the  inner  and  outer  layers  are  patched.  The  outer  layer  is  solved  in  a 
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truncated  computational  domain  with  the  initial  value  line  positioned 
slightly  upstream  of  the  start  of  the  interaction  region.  Thus,  the 
outer  solution  depends  on  two  geometric  parameters,  the  patch  distance 
and  the  axial  location  of  the  initial  value  line.  Numerical  tests  have 
shown  that  the  solution  is  relatively  insensitive  to  the  location  of  the 
initial  value  line  but  is  affected  significantly  by  the  patch  distance. 
Initial  calculations  indicate  that  the  patch  distance,  which  serves  as  a 
calibration  parameter,  depends  oniy  on  Reynolds  number.  Once  this 
parameter  was  determined,  computed  body  pressure  distributions  in  the 
viscid- inviscid  interaction  region  (for  three  test  bodies)  were  found  to 
follow  the  experimental  values  very  well. 
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"PO 


F 

G 

H 


static  pressure  coefficient 
total  head 

reduced  vorticity  -  defined  by  Eq.  (3) 

perturbation  stream  function  --  defined  by  Eq.  (2) 

derivative  of  G  with  respect  to  n  *  see  Eq.  (40) 

.....  .  -  oG 

spline  derivative  approximation  or  — 


,H 


(NO 


b  JE, 


Nn 

r 

rb 

Re 

u 


spline  derivative  approximation  of 

number  of  nodes  in  the  E, -direction 
number  of  5  nodes  on  the  body 
number  of  nodes  in  the  redirection 
radial  coordinate 
body  radius 

Reynolds  number,  based  on  body  length  and  free-stream  speed 

velocity  component  in  x-direction 

velocity  component  in  r-direction 

axial  coordinate 

vorticity  magnitude 

transformed  coordinate  along  body  and  centerline 
transformed  coordinate  away  from  body  and  centerline 
stream  function 

angle  of  tangent  to  body  surface 
meridian  angle 


All  other  quantities  are  defined  in  the  text. 
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Figure  Captions 

1.  C-Grip  Mapping. 

2.  Coordinates  in  Nose  Region. 

3.  Reduced  Vorticity  Distribution  on  Initial  Value  Line,  NSRDC  3odv. 

4.  Viscid-Inviscid  Interaction  Body  Pressure  Distribution,  F-57  Body. 

5.  Viscid-Inviscid  Interaction  Body  Pressure  Distribution,  Modified 
Spheroid. 

6.  Viscid-Inviscid  Interaction  Body  Pressure  Distribution ,  NSRLC  Body 
No.  1. 

7.  Comparison  of  Turbulent  Boundary-Layer  and  Frozen  Vorticity 
Velocity  Profiles  at  x  =  0.844,  F-57  Body. 

8.  Comparison  of  Turbulent  Boundary- Layer  and  Frozen  Vorticity 
Velocity  Profiles  at  x  =  0.933,  F-57  Body. 

9.  Comparison  of  Turbulent  Boundary-Layer  and  Frozen  Vorticity 
Velocity  Profiles  at  x  =  0.978,  F-57  Body. 
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